Case study for Covid-19 B-cells: downstream comparison and visualization

1 Load libraries

library(SingleCellExperiment)
library(scuttle)
library(ggplot2)
library(patchwork)
library(scater)
library(scran)
library(edgeR)
library(limma)
library(pbapply)
library(iCOBRA)
library(harmonicmeanp)
library(stageR)
library(dplyr)
library(ComplexHeatmap)

here::i_am("4_CovidCase_downstream.Rmd")

## here() starts at /Users/jg/Desktop/PhD/DD_project/DD_cases/Covid

2 Load DGE analysis data and results

pb_ct <- readRDS("./objects/pb_ct_filt.rds")
DGE_res <- readRDS("./objects/DGE_res_edgeR_NB.rds")

3 Load DD analysis data and results

pb_bin_ct <- readRDS("./objects/pb_bin_ct_filt.rds")
DD_res <- readRDS("./objects/DD_res_edgeR_NB_optim.rds")

4 Venn diagram

pdf("./figures/Venndiagram_DD_DE.pdf",
            width     = 5,
            height    = 3,
            pointsize = 4)
for (j in seq_along(DD_res)) {
  
  levels_status <- levels(pb_bin_ct$`B cell`$Status_on_day_collection_summary)
  
  par(mfrow=c(2,3))
  for(i in 1:5){
    pval_df <- data.frame(DGE = DGE_res[[j]][[i]]$PValue,
                          DD = DD_res[[j]][[i]]$PValue)
    rownames(pval_df) <- rownames(DGE_res[[j]][[i]])
    pad <- pval_df
    pad$DGE <- p.adjust(pad$DGE, method = "BH")
    pad$DD <- p.adjust(pad$DD, method = "BH")
    cobradata <- COBRAData(pval = pval_df, padj=pad)
    cobraperf <- calculate_performance(cobradata, splv = "none",maxsplit = 4)
    cobraplot <- prepare_data_for_plot(cobraperf, colorscheme = "Dark2", facetted = TRUE)
    
    plot_overlap(cobraplot, main = paste0(names(DD_res)[j], "\n", levels_status[i+1]),
                 mar = rep(2,4))
  }
}
dev.off()

## quartz_off_screen 
##                 2

5 Stagewise testing

stageR_res <- vector(mode="list",length = 5)
names(stageR_res) <- names(DD_res)

for (j in seq_along(DD_res)) {
  stageR_res_j <- vector(mode="list",length = 5)
  names(stageR_res_j) <- levels(pb_bin_ct$`B cell`$Status_on_day_collection_summary)[2:6]
  
  for(i in 1:5){
    pScreen <- rep(NA, nrow(DD_res[[j]][[i]]))
    for(h in 1:length(pScreen)){
        if(is.na(DD_res[[j]][[i]]$PValue[h])){
          pScreen[h] <- DGE_res[[j]][[i]]$PValue[h]
        } else {
          pScreen[h] <- hmp.stat(c(DD_res[[j]][[i]]$PValue[h],
                                   DGE_res[[j]][[i]]$PValue[h]),
                                 w=NULL)
        }
    }
    names(pScreen) <- rownames(DD_res[[j]][[i]])
    
    #confirmation stage
    pConfirmation <- as.matrix(cbind(DD_res[[j]][[i]]$PValue, DGE_res[[j]][[i]]$PValue))
    dimnames(pConfirmation) <- list(rownames(DD_res[[j]][[i]]), c("DD","DE"))
    
    # stageWise analysis
    stageRObj <- stageR(pScreen=pScreen, 
                        pConfirmation=pConfirmation, 
                        pScreenAdjusted=FALSE)
    stageRObj <- stageWiseAdjustment(object=stageRObj, 
                                     method="none", 
                                     alpha=0.05,
                                     allowNA = TRUE)
    res <- getResults(stageRObj)
    print(colSums(res)) #stage-wise analysis results
    
    if(colSums(res)[1] == 1){
      DD_DE_Sig_j_i <- getAdjustedPValues(stageRObj, 
                                      onlySignificantGenes=TRUE, 
                                      order=FALSE)
    } else{
      DD_DE_Sig_j_i <- getAdjustedPValues(stageRObj, 
                                      onlySignificantGenes=TRUE, 
                                      order=TRUE)
    }
    
    if(is.null(DD_DE_Sig_j_i)){
      DD_DE_Sig_j_i <- NA
    }
    
    stageR_res_j[[i]] <- DD_DE_Sig_j_i
  }
  stageR_res[[j]] <- stageR_res_j
}

## padjScreen         DD         DE 
##          0          0          0 
## padjScreen         DD         DE 
##          1          1          0 
## padjScreen         DD         DE 
##          5          3          5 
## padjScreen         DD         DE 
##          1          1          1 
## padjScreen         DD         DE 
##          0          0          0 
## padjScreen         DD         DE 
##          0          0          0 
## padjScreen         DD         DE 
##         95         33         81 
## padjScreen         DD         DE 
##        384        213        326 
## padjScreen         DD         DE 
##        180         72        138 
## padjScreen         DD         DE 
##         15         10          6 
## padjScreen         DD         DE 
##          1          1          1 
## padjScreen         DD         DE 
##         71         40         55 
## padjScreen         DD         DE 
##        326        154        297 
## padjScreen         DD         DE 
##         88         52         74 
## padjScreen         DD         DE 
##         23         19          9 
## padjScreen         DD         DE 
##          0          0          0 
## padjScreen         DD         DE 
##        266        161        229 
## padjScreen         DD         DE 
##       2345       1580       2239 
## padjScreen         DD         DE 
##       1529        812       1474 
## padjScreen         DD         DE 
##       1149        453       1125 
## padjScreen         DD         DE 
##          0          0          0 
## padjScreen         DD         DE 
##         36         21         21 
## padjScreen         DD         DE 
##        232        162        159 
## padjScreen         DD         DE 
##         99         30         84 
## padjScreen         DD         DE 
##         16         14          3

for(j in 1:5){
  for(i in 1:5){
    if(all(is.na(stageR_res[[j]][[i]]))){
      stagewise <- 0
    } else if(length(stageR_res[[j]][[i]])==3){
      stagewise <- 1
    } else {
      stagewise <- sum(stageR_res[[j]][[i]][,1] < 0.05, na.rm=TRUE)
    }
    sig_DD <- rownames(DD_res[[j]][[i]])[which(p.adjust(DD_res[[j]][[i]]$PValue, method="BH") < 0.05)]
    sig_DE <- rownames(DGE_res[[j]][[i]])[which(p.adjust(DGE_res[[j]][[i]]$PValue, method="BH") < 0.05)]

    separate <- length(unique(sig_DD, sig_DE))
    
    print(c(stagewise, separate))
  }
}

## [1] 0 0
## [1] 1 2
## [1] 5 3
## [1] 1 1
## [1] 0 6
## [1] 0 0
## [1] 95  8
## [1] 384 191
## [1] 180  50
## [1] 15 11
## [1] 1 1
## [1] 71 36
## [1] 326 110
## [1] 88 43
## [1] 23 27
## [1] 0 0
## [1] 266 135
## [1] 2345 1417
## [1] 1529  555
## [1] 1149  224
## [1] 0 0
## [1] 36 19
## [1] 232 164
## [1] 99 11
## [1] 16 50

6 Visualize offset

par(mfrow=c(2,3))
for (element in pb_bin_ct) {
    bin_counts <- assay(element)

    of <- colMeans(sweep(bin_counts, 2, element$ncells, "/")) 
    logitOf <- log(of/(1-of))
    
    plot(x = jitter(as.numeric(element$Status_on_day_collection_summary)), 
         y = of,
         col = element$Site, 
         pch = 19, 
         cex = 0.8, 
         main = "Offset per Status, colored by Site",
         xlab = "Covid status")
}

7 Visualization helper functions

visualize_DGE <- function(pb, genes){
    cd <- colData(pb)
    cd <- droplevels(cd)
    libSize <- pb %>%
        counts %>%
        colSums %>%
        unname

    gg_list <- lapply(genes, function(gene){
      
        data <- data.frame(cpm = log2(counts(pb)[gene,]+0.5) - log2(libSize+1)+log2(1e6),
                           status = cd$Status_on_day_collection_summary,
                           batch = paste0(cd$Site,"_",cd$sex))
        data$cpm[is.infinite(data$cpm)] <- NA
        
        levels(data$batch) <- c("Cambridge_female", "Cambridge_male", 
                                "Ncl_female", "Ncl_male")
        
        data$status_batch <- factor(paste0(data$status, "_", data$batch),
                                    levels = paste0(rep(levels(data$status),each=4),
                                                    "_", 
                                                    rep(levels(data$batch), times=2)))
        data$dotposition <- data$status_batch
        levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
                                         "0.9" = levels(data$status_batch)[2],
                                         "1.1" = levels(data$status_batch)[3],
                                         "1.3" = levels(data$status_batch)[4],
                                         "1.7" = levels(data$status_batch)[5],
                                         "1.9" = levels(data$status_batch)[6],
                                         "2.1" = levels(data$status_batch)[7],
                                         "2.3" = levels(data$status_batch)[8])
        data$dotposition <- as.numeric(as.character(data$dotposition))
      
        gg <- ggplot(data = data, aes(x=status,y=cpm)) +
            geom_violin() +
            geom_jitter(aes(x=dotposition, 
                          y=cpm, 
                          col=batch),
                      size=2, 
                      width = 0.05) +
            ggtitle(paste0("Log CPM: ",gene)) +
            theme_bw() +
            theme(plot.title = element_text(size=11))
        return(gg)
    })
    return(gg_list)
}

visualize_DD_LOR <- function(pb, genes){
    bin_counts <- assay(pb)
    cd <- colData(pb)
    cd <- droplevels(cd)
    of <- colMeans(sweep(bin_counts, 2, cd$ncells, "/")) 
    odds_global <- of/(1-of)

    gg_list <- lapply(genes, function(gene){
      odds_gene <- bin_counts[gene,]/(cd$ncells - bin_counts[gene,])
      
      data <- data.frame(odds_gene = odds_gene,
                         odds_global = odds_global,
                         LOR = log(odds_gene/odds_global),
                         status = cd$Status_on_day_collection_summary,
                         batch = paste0(cd$Site,"_",cd$sex))
      data$LOR[is.infinite(data$LOR)] <- NA
      
      levels(data$batch) <- c("Cambridge_female", "Cambridge_male", 
                                "Ncl_female", "Ncl_male")
        
      data$status_batch <- factor(paste0(data$status, "_", data$batch),
                                  levels = paste0(rep(levels(data$status),each=4),
                                                  "_", 
                                                  rep(levels(data$batch), times=2)))
      data$dotposition <- data$status_batch
      levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
                                       "0.9" = levels(data$status_batch)[2],
                                       "1.1" = levels(data$status_batch)[3],
                                       "1.3" = levels(data$status_batch)[4],
                                       "1.7" = levels(data$status_batch)[5],
                                       "1.9" = levels(data$status_batch)[6],
                                       "2.1" = levels(data$status_batch)[7],
                                       "2.3" = levels(data$status_batch)[8])
      data$dotposition <- as.numeric(as.character(data$dotposition))
      
      gg <- ggplot(data = data,aes(x=status, y= LOR)) +
          geom_violin() +
          geom_jitter(aes(x=dotposition, 
                          y=LOR, 
                          col=batch),
                      size=2, 
                      width = 0.05) +
          ggtitle(paste0("LOR: ", gene)) +
          theme_bw() +
          theme(plot.title = element_text(size=11))
      return(gg)
    })
    return(gg_list)
}

visualize_DD_logprop <- function(pb, genes){
    bin_counts <- assay(pb)
    cd <- colData(pb)
    cd <- droplevels(cd)
    logprop <- log(sweep(bin_counts, 2, cd$ncells, "/"))

    gg_list <- lapply(genes, function(gene){

      data <- data.frame(logprop = logprop[gene,],
                         status = cd$Status_on_day_collection_summary,
                         batch = paste0(cd$Site,"_",cd$sex))
      data$logprop[is.infinite(data$logprop)] <- NA
      
      levels(data$batch) <- c("Cambridge_female", "Cambridge_male", 
                                "Ncl_female", "Ncl_male")
        
      data$status_batch <- factor(paste0(data$status, "_", data$batch),
                                  levels = paste0(rep(levels(data$status),each=4),
                                                  "_", 
                                                  rep(levels(data$batch), times=2)))
      data$dotposition <- data$status_batch
      levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
                                       "0.9" = levels(data$status_batch)[2],
                                       "1.1" = levels(data$status_batch)[3],
                                       "1.3" = levels(data$status_batch)[4],
                                       "1.7" = levels(data$status_batch)[5],
                                       "1.9" = levels(data$status_batch)[6],
                                       "2.1" = levels(data$status_batch)[7],
                                       "2.3" = levels(data$status_batch)[8])
      data$dotposition <- as.numeric(as.character(data$dotposition))
      
      gg <- ggplot(data = data,aes(x=status, y= logprop)) +
          geom_violin() +
          geom_jitter(aes(x=dotposition, 
                          y=logprop, 
                          col=batch),
                      size=2, 
                      width = 0.05) +
          ggtitle(paste0("logprop: ", gene)) +
          theme_bw() +
          theme(plot.title = element_text(size=11))
      return(gg)
    })
    return(gg_list)
}

visualize_DD_prop <- function(pb, genes){
    bin_counts <- assay(pb)
    cd <- colData(pb)
    cd <- droplevels(cd)
    proportions <- sweep(bin_counts, 2, cd$ncells, "/")

    gg_list <- lapply(genes, function(gene){

        data <- data.frame(proportion = proportions[gene,],
                           status = cd$Status_on_day_collection_summary,
                           batch = paste0(cd$Site,"_",cd$sex))
        data$proportion[is.infinite(data$proportion)] <- NA
        
        levels(data$batch) <- c("Cambridge_female", "Cambridge_male", 
                                "Ncl_female", "Ncl_male")
        
        data$status_batch <- factor(paste0(data$status, "_", data$batch),
                                    levels = paste0(rep(levels(data$status),each=4),
                                             "_", 
                                             rep(levels(data$batch), times=2)))
        
        data$dotposition <- data$status_batch
        levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
                                         "0.9" = levels(data$status_batch)[2],
                                         "1.1" = levels(data$status_batch)[3],
                                         "1.3" = levels(data$status_batch)[4],
                                         "1.7" = levels(data$status_batch)[5],
                                         "1.9" = levels(data$status_batch)[6],
                                         "2.1" = levels(data$status_batch)[7],
                                         "2.3" = levels(data$status_batch)[8])
        data$dotposition <- as.numeric(as.character(data$dotposition))
        
        gg <- ggplot(data = data,aes(x=status, y= proportion)) +
            geom_violin() +
            geom_jitter(aes(x=dotposition, 
                            y=proportion, 
                            col=batch),
                        size=2, 
                        width = 0.05) +
            ylim(-0.01,1.01) + 
            ggtitle(paste0("prop: ", gene)) +
            theme_bw() +
            theme(plot.title = element_text(size=11))
      return(gg)
    })
    return(gg_list)
}

7.1 Prepare data for plotting

pb_counts_list <- list(Naive_moderate = pb_ct$`naive B cell`[,which(pb_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
                       Naive_critical = pb_ct$`naive B cell`[,which(pb_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))],
                       Unswitched_moderate = pb_ct$`unswitched memory B cell`[,which(pb_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
                       Unswitched_critical = pb_ct$`unswitched memory B cell`[,which(pb_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))])

pb_bin_list <- list(Naive_moderate = pb_bin_ct$`naive B cell`[,which(pb_bin_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
                    Naive_critical = pb_bin_ct$`naive B cell`[,which(pb_bin_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))],
                    Unswitched_moderate = pb_bin_ct$`unswitched memory B cell`[,which(pb_bin_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
                    Unswitched_critical = pb_bin_ct$`unswitched memory B cell`[,which(pb_bin_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))])

7.2 Naive B cells - moderate

7.2.1 Top 10 DD

##              logFC   logCPM        F       PValue
## TRBC2   -0.4761722 8.770894 62.98263 1.584941e-12
## ICAM3   -0.5147748 8.565486 60.72267 3.335492e-12
## TGFBR2  -0.5887690 7.925526 59.98712 4.258692e-12
## RNF141  -0.6553944 7.463683 59.21566 5.509029e-12
## MKI67    4.0945860 4.318984 59.42586 5.640242e-12
## SLC43A2 -0.9289432 5.793013 58.13510 7.916771e-12
## MYLIP   -0.7271692 7.722048 55.98090 2.607019e-11
## ADD1    -0.4508512 7.924785 52.87000 4.796916e-11
## BLCAP   -0.6159960 7.496801 51.76362 7.057862e-11
## RRM2     4.6558042 4.262710 51.41045 1.163980e-10

7.2.2 Top 10 DD not DGE (pvalue)

##               logFC   logCPM        F       PValue
## ATP5F1A  -0.3454894 8.359434 28.74228 4.368316e-07
## IGHV6-1   1.1558553 4.874379 18.74248 3.229214e-05
## ACTG1    -0.3194227 9.447387 18.32626 4.351886e-05
## GDI2     -0.2332617 8.832662 17.94472 4.629560e-05
## HSPA8    -0.2259126 9.136426 15.85318 1.208046e-04
## HMGB2    -0.3538589 7.858732 14.95412 1.925820e-04
## LITAF    -0.2540148 8.096516 14.50821 2.264128e-04
## IGLV1-51  0.9283355 6.242311 14.45082 2.588274e-04
## SHMT2    -0.2586742 8.166360 13.93754 2.964212e-04
## PKM      -0.2190923 8.514743 12.81372 5.065802e-04

7.2.3 Manuscript figure

pdf("./figures/figure3.pdf",
            width     = 8,
            height    = 5,
            pointsize = 4)
DD_prop_plots[[3]] + DGE_plots[[3]] + plot_layout(guides = "collect")
dev.off()

## quartz_off_screen 
##                 2

7.2.4 Top 10 DD not DGE (lfc)

##                logFC   logCPM         F       PValue
## LINC01238 -1.2922121 2.983519  7.841155 5.999375e-03
## SYNC      -0.4995439 5.612165  7.935728 5.713336e-03
## PIK3R5    -0.4212666 5.776195  9.260385 2.906765e-03
## HNRNPUL2  -0.3752241 5.923456  7.463781 7.297005e-03
## HMGB2     -0.3538589 7.858732 14.954121 1.925820e-04
## ATP5F1A   -0.3454894 8.359434 28.742279 4.368316e-07
## SAMM50    -0.3201015 6.653094  8.159834 5.090509e-03
## ACTG1     -0.3194227 9.447387 18.326257 4.351886e-05
## DDT       -0.2913479 8.224187 10.601145 1.551812e-03
## UBA1      -0.2710409 7.063174  7.961008 5.639294e-03

7.3 Naive B cells - critical

7.3.1 Top 10 DD

##             logFC   logCPM        F       PValue
## MKI67   4.0121620 4.318984 44.40367 1.050734e-09
## CDKN2D -0.6837345 7.275963 33.51255 6.339351e-08
## PSAT1   3.3397088 3.811498 31.77674 1.269090e-07
## TRBC2  -0.4243709 8.770894 31.51889 1.408042e-07
## RRM2    4.2225069 4.262710 30.89135 2.146758e-07
## BIRC5   3.8207565 4.105887 30.52895 2.319837e-07
## UBE2C   4.3883574 3.738962 27.41210 7.581503e-07
## SEMA4B -0.6655763 7.226843 26.10699 1.367672e-06
## PCLAF   2.9412461 4.277981 25.87200 1.562934e-06
## BUB1    2.6976148 4.102725 25.68458 1.565488e-06

7.3.2 Top 10 DD not DGE (pvalue)

##               logFC   logCPM        F       PValue
## IFITM1    0.7048667 9.359272 15.57898 0.0001553255
## IGHV6-1   1.1491385 4.874379 13.20591 0.0004198879
## MCM5     -0.3864310 7.280573 12.76297 0.0005190805
## SNX6      0.2209095 8.361851 12.55092 0.0005748248
## MYC       0.4780873 8.027438 11.92146 0.0008381220
## S100A8    1.1834574 8.245194 11.59541 0.0009842203
## IGKV1D-8  1.1762445 4.358691 11.28444 0.0010632792

7.3.3 Top 10 DD not DGE (lfc)

##               logFC   logCPM        F       PValue
## MCM5     -0.3864310 7.280573 12.76297 0.0005190805
## SNX6      0.2209095 8.361851 12.55092 0.0005748248
## MYC       0.4780873 8.027438 11.92146 0.0008381220
## IFITM1    0.7048667 9.359272 15.57898 0.0001553255
## IGHV6-1   1.1491385 4.874379 13.20591 0.0004198879
## IGKV1D-8  1.1762445 4.358691 11.28444 0.0010632792
## S100A8    1.1834574 8.245194 11.59541 0.0009842203

7.4 Unswitched memory B cells - moderate

7.4.1 Top 10 DD not DGE (pvalue)

##               logFC   logCPM        F       PValue
## PSME2     0.5507604 9.476988 23.54312 5.383867e-06
## IGKC      0.6453526 9.328480 22.97122 9.375926e-06
## PMAIP1    1.1439169 8.479374 21.65406 2.421569e-05
## KIAA0040 -0.7636898 8.622081 19.16872 3.364242e-05
## NFKBID    1.0154309 8.356222 19.95561 3.416495e-05
## DDX21     0.4502855 9.322344 17.97053 5.154182e-05
## EIF2AK2   1.2710702 8.693864 19.70351 7.317991e-05
## CD27     -0.7000375 8.614283 17.48454 7.418096e-05
## SMIM14   -0.4986694 9.076394 16.29166 1.090757e-04
## ARHGAP30 -0.5941593 8.723185 16.06569 1.207913e-04

7.4.2 Top 10 DD not DGE (lfc)

##               logFC   logCPM         F       PValue
## GPM6A    -1.5599708 7.348838 13.811090 0.0004432711
## C3orf38  -1.1731767 7.454788 12.955172 0.0005066061
## MID1IP1  -1.1504175 7.432318 10.821879 0.0014068824
## OSBPL10  -0.9746514 7.706646 13.537482 0.0003856664
## BCL7A    -0.8985348 7.871920 12.755253 0.0005953312
## TP53I11  -0.8918054 7.622110  9.614723 0.0025329619
## GPAA1    -0.8884655 7.861906 14.407231 0.0002576068
## RNF113A  -0.8447010 7.695519  9.766845 0.0023496663
## TSPYL1   -0.8150866 8.014314 14.454007 0.0002521083
## ALDH16A1 -0.8027842 7.974603 11.063781 0.0013002066

7.5 Unswitched memory B cells - critical

7.5.1 Top 10 DD not DGE (pvalue)

##             logFC    logCPM        F       PValue
## EIF5A   0.6137839  9.754504 32.12061 1.512792e-07
## DDX21   0.6480403  9.322344 27.12606 1.083771e-06
## IGKC    0.8428443  9.328480 28.18164 1.322635e-06
## COX6A1  0.4530417  9.825657 21.02758 1.363526e-05
## CD27   -0.9792277  8.614283 19.43624 3.235216e-05
## RHOA    0.4199441 10.007795 18.86143 3.482467e-05
## CD1C   -1.2118898  8.643620 20.93737 3.570295e-05
## ZNF581 -1.1502640  8.320019 18.66492 3.795786e-05
## TNIP2   1.2475800  7.881964 18.31475 4.427646e-05
## PSME2   0.5870198  9.476988 18.50603 4.451462e-05

7.5.2 Top 10 DD not DGE (lfc)

##                logFC   logCPM        F       PValue
## LCN8      -5.9292313 6.824761 13.91315 6.293238e-04
## GPM6A     -2.2523892 7.348838 13.05752 6.154990e-04
## FCGRT     -1.2570043 7.717782 13.08554 4.765082e-04
## CD1C      -1.2118898 8.643620 20.93737 3.570295e-05
## ZNF581    -1.1502640 8.320019 18.66492 3.795786e-05
## CD27      -0.9792277 8.614283 19.43624 3.235216e-05
## RIN3      -0.7505344 8.341573 12.68901 5.742842e-04
## FXR1      -0.7379673 8.517854 12.39759 6.591417e-04
## FCRL2     -0.7276242 8.848491 15.26450 1.795342e-04
## LINC01857 -0.6606280 9.070504 13.87236 3.627673e-04

8 Single-cell plots

sce <- readRDS("./objects/sce_Covid_Bcells.rds")

sce_sub <- sce[,sce$cell_type_curated == "naive B cell"]
sce_sub <- sce_sub[,sce_sub$donor_id %in% pb_bin_list$Naive_moderate$donor_id]

lnc <- logNormCounts(sce_sub,name="logNormCounts")
lnc$Site_sex <- paste0(lnc$Site, "_", lnc$sex)

gene <- "ACTG1"
DD_prop_plots <- visualize_DD_prop(pb = pb_bin_list$Naive_moderate, 
                                   genes = gene)
DGE_plots <- visualize_DGE(pb = pb_counts_list$Naive_moderate, 
                           genes = gene)

print(DD_prop_plots[[1]] + DGE_plots[[1]] + plot_layout(guides = "collect"))

df_hlp <- data.frame(counts = assays(sce_sub)$counts[gene,],
                     lnc = assays(lnc)$logNormCounts[gene,],
                     Site_sex = lnc$Site_sex,
                     Status = lnc$Status_on_day_collection_summary,
                     Donor = as.factor(as.character(lnc$donor_id)))

df_hlp <- df_hlp[df_hlp$Site_sex == "Ncl_male",]

gg1 <- ggplot(data = df_hlp, aes(x=Donor,y=counts, fill=Status)) +
            geom_violin() +
            facet_grid(~Status, space = "free_x", scales = "free_x") +
            ggtitle(paste0("Counts: ",gene)) +
            theme_bw() +
            theme(plot.title = element_text(size=11),
                  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

gg2 <- ggplot(data = df_hlp, aes(x=Donor,y=lnc, fill=Status)) +
            geom_violin() +
            facet_grid(~Status, space = "free_x", scales = "free_x") +
            ggtitle(paste0("LogNormCounts: ",gene)) +
            theme_bw() +
            theme(plot.title = element_text(size=11),
                  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

pdf("./figures/figure4.pdf",
            width     = 8,
            height    = 5,
            pointsize = 4)
gg2
dev.off()

## quartz_off_screen 
##                 2

9 Get percentage 0-1-2 for single-cell level data

pb_ct <- readRDS("./objects/pb_ct_filt.rds")
sce <- readRDS("./objects/sce_Covid_Bcells.rds")

# B cell
sce_hlp <- sce[,sce$cell_type_curated == "B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`B cell`$donor_id]
a_entries <- nrow(sce_hlp)*ncol(sce_hlp)
a_0 <- sum(assay(sce_hlp) == 0)
a_1 <- sum(assay(sce_hlp) == 1)
a_2 <- sum(assay(sce_hlp) == 2)
a_r <- sum(assay(sce_hlp) > 2)

# class switched memory B cell
sce_hlp <- sce[,sce$cell_type_curated == "class switched memory B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`class switched memory B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`class switched memory B cell`$donor_id]
b_entries <- nrow(sce_hlp)*ncol(sce_hlp)
b_0 <- sum(assay(sce_hlp) == 0)
b_1 <- sum(assay(sce_hlp) == 1)
b_2 <- sum(assay(sce_hlp) == 2)
b_r <- sum(assay(sce_hlp) > 2)

# immature B cell
sce_hlp <- sce[,sce$cell_type_curated == "immature B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`immature B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`immature B cell`$donor_id]
c_entries <- nrow(sce_hlp)*ncol(sce_hlp)
c_0 <- sum(assay(sce_hlp) == 0)
c_1 <- sum(assay(sce_hlp) == 1)
c_2 <- sum(assay(sce_hlp) == 2)
c_r <- sum(assay(sce_hlp) > 2)

# naive B cell
sce_hlp <- sce[,sce$cell_type_curated == "naive B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`naive B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`naive B cell`$donor_id]
d_entries <- nrow(sce_hlp)*ncol(sce_hlp)
d_0 <- sum(assay(sce_hlp) == 0)
d_1 <- sum(assay(sce_hlp) == 1)
d_2 <- sum(assay(sce_hlp) == 2)
d_r <- sum(assay(sce_hlp) > 2)

# unswitched memory B cell
sce_hlp <- sce[,sce$cell_type_curated == "unswitched memory B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`unswitched memory B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`unswitched memory B cell`$donor_id]
e_entries <- nrow(sce_hlp)*ncol(sce_hlp)
e_0 <- sum(assay(sce_hlp) == 0)
e_1 <- sum(assay(sce_hlp) == 1)
e_2 <- sum(assay(sce_hlp) == 2)
e_r <- sum(assay(sce_hlp) > 2)

(a_0+b_0+c_0+d_0+e_0)/(a_entries+b_entries+c_entries+d_entries+e_entries)*100

## [1] 86.09558

(a_1+b_1+c_1+d_1+e_1)/(a_entries+b_entries+c_entries+d_entries+e_entries)*100

## [1] 9.386102

(a_2+b_2+c_2+d_2+e_2)/(a_entries+b_entries+c_entries+d_entries+e_entries)*100

## [1] 1.924958

(a_r+b_r+c_r+d_r+e_r)/(a_entries+b_entries+c_entries+d_entries+e_entries)*100

## [1] 2.593364

10 Session info

sessionInfo()

## R version 4.2.0 alpha (2022-04-04 r82084)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ComplexHeatmap_2.14.0       dplyr_1.1.2                
##  [3] stageR_1.20.0               harmonicmeanp_3.0          
##  [5] FMStable_0.1-4              iCOBRA_1.26.0              
##  [7] pbapply_1.7-2               edgeR_3.40.2               
##  [9] limma_3.54.2                scran_1.26.2               
## [11] scater_1.26.1               patchwork_1.1.2            
## [13] ggplot2_3.4.2               scuttle_1.8.4              
## [15] SingleCellExperiment_1.20.1 SummarizedExperiment_1.28.0
## [17] Biobase_2.58.0              GenomicRanges_1.50.2       
## [19] GenomeInfoDb_1.34.9         IRanges_2.32.0             
## [21] S4Vectors_0.36.2            BiocGenerics_0.44.0        
## [23] MatrixGenerics_1.10.0       matrixStats_1.0.0          
## 
## loaded via a namespace (and not attached):
##   [1] ggbeeswarm_0.7.2          colorspace_2.1-0         
##   [3] rjson_0.2.21              ellipsis_0.3.2           
##   [5] rprojroot_2.0.3           circlize_0.4.15          
##   [7] bluster_1.8.0             XVector_0.38.0           
##   [9] GlobalOptions_0.1.2       BiocNeighbors_1.16.0     
##  [11] clue_0.3-64               rstudioapi_0.14          
##  [13] farver_2.1.1              ggrepel_0.9.3            
##  [15] DT_0.28                   fansi_1.0.4              
##  [17] codetools_0.2-19          sparseMatrixStats_1.10.0 
##  [19] doParallel_1.0.17         cachem_1.0.8             
##  [21] knitr_1.43                jsonlite_1.8.7           
##  [23] cluster_2.1.4             png_0.1-8                
##  [25] shinydashboard_0.7.2      shiny_1.7.4.1            
##  [27] compiler_4.2.0            dqrng_0.3.0              
##  [29] Matrix_1.5-4.1            fastmap_1.1.1            
##  [31] cli_3.6.1                 later_1.3.1              
##  [33] BiocSingular_1.14.0       htmltools_0.5.5          
##  [35] tools_4.2.0               rsvd_1.0.5               
##  [37] igraph_1.5.0.1            gtable_0.3.3             
##  [39] glue_1.6.2                GenomeInfoDbData_1.2.9   
##  [41] reshape2_1.4.4            Rcpp_1.0.11              
##  [43] jquerylib_0.1.4           vctrs_0.6.3              
##  [45] iterators_1.0.14          DelayedMatrixStats_1.20.0
##  [47] xfun_0.39                 stringr_1.5.0            
##  [49] beachmat_2.14.2           mime_0.12                
##  [51] lifecycle_1.0.3           irlba_2.3.5.1            
##  [53] statmod_1.5.0             zlibbioc_1.44.0          
##  [55] scales_1.2.1              promises_1.2.0.1         
##  [57] shinyBS_0.61.1            parallel_4.2.0           
##  [59] RColorBrewer_1.1-3        yaml_2.3.7               
##  [61] gridExtra_2.3             UpSetR_1.4.0             
##  [63] sass_0.4.7                stringi_1.7.12           
##  [65] highr_0.10                foreach_1.5.2            
##  [67] ScaledMatrix_1.6.0        BiocParallel_1.32.6      
##  [69] shape_1.4.6               rlang_1.1.1              
##  [71] pkgconfig_2.0.3           bitops_1.0-7             
##  [73] evaluate_0.21             lattice_0.21-8           
##  [75] ROCR_1.0-11               labeling_0.4.2           
##  [77] htmlwidgets_1.6.2         tidyselect_1.2.0         
##  [79] here_1.0.1                plyr_1.8.8               
##  [81] magrittr_2.0.3            R6_2.5.1                 
##  [83] generics_0.1.3            metapod_1.6.0            
##  [85] DelayedArray_0.24.0       DBI_1.1.3                
##  [87] pillar_1.9.0              withr_2.5.0              
##  [89] RCurl_1.98-1.12           tibble_3.2.1             
##  [91] crayon_1.5.2              utf8_1.2.3               
##  [93] rmarkdown_2.23            viridis_0.6.4            
##  [95] GetoptLong_1.0.5          locfit_1.5-9.8           
##  [97] digest_0.6.33             xtable_1.8-4             
##  [99] httpuv_1.6.11             munsell_0.5.0            
## [101] beeswarm_0.4.0            viridisLite_0.4.2        
## [103] vipor_0.4.5               bslib_0.5.0